stanEx4
execute mcmc sampling
goMCMC=function(mdl,data,smp=500,wrm=100,th=1){
mcmc=mdl$sample(
data=data,
seed=1,
chains=4,
iter_sampling=smp,
iter_warmup=wrm,
thin=th,
refresh=1000
)
mcmc
}see mcmc result and parameters
seeMCMC=function(mcmc,exc='',ch=T){ # exclude 'exc' parameters from seeing
print(mcmc)
prs=mcmc$metadata()$model_params[-1] # reject lp__
for(pr in prs){
if(grepl('^y',pr)) next # not show predictive value "y*" information
if(exc!='' && grepl(paste0('^',exc),pr)) next
drw=mcmc$draws(pr)
if(ch){
par(mfrow=c(1,4))
drw[,1,] |> plot(type='l',main='chain1',ylab=pr)
drw[,2,] |> plot(type='l',main='chain2',ylab=pr)
drw[,3,] |> plot(type='l',main='chain3',ylab=pr)
drw[,4,] |> plot(type='l',main='chain4',ylab=pr)
}
par(mfrow=c(1,2))
drw |> hist(main=pr,xlab='')
drw |> density() |> plot(main=pr)
}
}fn=function(mdl,data,smp=500,wrm=100){
mle=mdl$optimize(data=data)
print(mle)
mcmc=goMCMC(mdl,data,smp,wrm)
mcmc$metadata()$stan_variables
seeMCMC(mcmc,'m')
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
}ex9
explanatory variable obs.x also has noise
x~N(x0.sx)
y~N(b0+b1*x0,s) ex8-0.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}ex9.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x10;
}
parameters{
real b0;
real b1;
real<lower=0> s;
real<lower=0> sx;
vector[N] x0;
}
model{
for(i in 1:N){
x[i]~normal(x0[i],sx);
y[i]~normal(b0+b1*x0[i],s);
}
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x0[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] x1;
vector[N1] y1;
for(i in 1:N1){
x1[i]=normal_rng(x10[i],sx);
m1[i]=b0+b1*x10[i];
y1[i]=normal_rng(m1[i],s);
}
}n=20
x0=runif(n,0,20)
x=rnorm(n,x0,2)
y=rnorm(n,x0*2+5,2)
qplot(x,y)
n1=10
#explanatory variable do not has noise
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -131976
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## 26 -37.5655 0.000105984 0.000189047 0.9582 0.9582 69
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.2 seconds.
## variable estimate
## lp__ -37.57
## b0 5.31
## b1 1.88
## s 3.97
## m0[1] 16.83
## m0[2] 14.97
## m0[3] 14.78
## m0[4] 11.97
## m0[5] 28.20
## m0[6] 27.81
## m0[7] 12.88
## m0[8] 24.81
## m0[9] 42.98
## m0[10] 5.24
## m0[11] 31.82
## m0[12] 20.31
## m0[13] 5.59
## m0[14] 37.79
## m0[15] 21.00
## m0[16] 18.53
## m0[17] 10.80
## m0[18] 14.35
## m0[19] 7.59
## m0[20] 37.71
## y0[1] 16.18
## y0[2] 8.69
## y0[3] 9.34
## y0[4] 13.05
## y0[5] 30.72
## y0[6] 29.23
## y0[7] 15.28
## y0[8] 24.01
## y0[9] 47.22
## y0[10] -1.53
## y0[11] 23.61
## y0[12] 25.20
## y0[13] 4.39
## y0[14] 32.02
## y0[15] 17.00
## y0[16] 25.43
## y0[17] 13.60
## y0[18] 6.00
## y0[19] 10.32
## y0[20] 36.68
## m1[1] 5.24
## m1[2] 9.43
## m1[3] 13.62
## m1[4] 17.82
## m1[5] 22.01
## m1[6] 26.21
## m1[7] 30.40
## m1[8] 34.59
## m1[9] 38.79
## m1[10] 42.98
## y1[1] 10.62
## y1[2] 8.09
## y1[3] 9.01
## y1[4] 19.50
## y1[5] 17.65
## y1[6] 25.88
## y1[7] 29.74
## y1[8] 40.45
## y1[9] 45.08
## y1[10] 42.03
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.4 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -37.75 -37.41 1.32 1.03 -40.29 -36.36 1.01 567 945
## b0 5.14 5.13 1.67 1.65 2.41 7.89 1.01 370 382
## b1 1.90 1.89 0.17 0.17 1.62 2.18 1.00 509 677
## s 4.49 4.37 0.83 0.75 3.41 5.93 1.00 1335 1120
## m0[1] 16.75 16.72 1.06 1.02 15.08 18.45 1.01 566 717
## m0[2] 14.87 14.85 1.11 1.08 13.08 16.64 1.01 479 624
## m0[3] 14.69 14.67 1.12 1.09 12.87 16.47 1.01 473 604
## m0[4] 11.85 11.83 1.25 1.22 9.80 13.85 1.01 409 435
## m0[5] 28.20 28.19 1.28 1.19 26.10 30.30 1.01 2125 1240
## m0[6] 27.80 27.78 1.26 1.19 25.73 29.89 1.01 2128 1240
## m0[7] 12.77 12.75 1.20 1.17 10.79 14.70 1.01 424 460
## m0[8] 24.79 24.77 1.12 1.06 23.02 26.68 1.01 1865 1214
## m0[9] 43.09 43.05 2.36 2.28 39.26 46.88 1.00 1054 1219
## m0[10] 5.07 5.05 1.68 1.66 2.32 7.83 1.01 370 382
## m0[11] 31.85 31.86 1.50 1.38 29.40 34.22 1.01 1861 1347
## m0[12] 20.25 20.22 1.02 1.00 18.63 21.98 1.01 932 1054
## m0[13] 5.43 5.41 1.65 1.63 2.72 8.12 1.01 371 385
## m0[14] 37.86 37.87 1.94 1.83 34.63 40.95 1.00 1351 1325
## m0[15] 20.95 20.91 1.03 1.01 19.32 22.68 1.01 1055 1094
## m0[16] 18.46 18.43 1.03 1.02 16.82 20.15 1.01 700 888
## m0[17] 10.67 10.66 1.31 1.28 8.51 12.76 1.01 395 390
## m0[18] 14.25 14.22 1.14 1.10 12.39 16.06 1.01 459 562
## m0[19] 7.44 7.41 1.51 1.51 4.97 9.87 1.01 375 386
## m0[20] 37.78 37.78 1.93 1.83 34.56 40.86 1.00 1357 1363
## y0[1] 16.91 16.85 4.82 4.52 8.99 24.79 1.00 1907 1631
## y0[2] 14.73 14.88 4.51 4.45 7.13 21.98 1.00 1480 1686
## y0[3] 14.60 14.57 4.81 4.55 6.58 22.62 1.00 1715 1790
## y0[4] 11.71 11.72 4.68 4.61 3.86 19.39 1.00 1626 1698
## y0[5] 28.27 28.39 4.61 4.30 20.77 35.79 1.00 2045 1977
## y0[6] 27.81 27.81 4.77 4.50 19.99 35.54 1.00 1878 1861
## y0[7] 12.72 12.63 4.69 4.46 5.14 20.24 1.00 1673 1784
## y0[8] 24.89 24.88 4.70 4.42 17.11 32.63 1.00 2030 2013
## y0[9] 42.96 42.85 5.28 5.12 34.64 51.87 1.01 1987 1832
## y0[10] 4.96 4.98 4.78 4.81 -2.97 12.65 1.00 1330 1671
## y0[11] 31.81 31.73 4.80 4.79 24.08 39.53 1.00 2042 1821
## y0[12] 20.29 20.21 4.66 4.45 12.69 27.82 1.00 2150 1864
## y0[13] 5.56 5.39 4.82 4.77 -2.34 13.42 1.00 1478 1805
## y0[14] 37.93 37.81 4.90 4.60 30.13 46.13 1.00 1997 1828
## y0[15] 20.94 20.99 4.76 4.56 13.25 28.68 1.00 1535 1522
## y0[16] 18.42 18.37 4.75 4.57 10.74 25.81 1.00 1734 1846
## y0[17] 10.61 10.49 4.73 4.47 3.28 18.39 1.00 1166 1363
## y0[18] 14.14 14.10 4.60 4.44 6.58 21.68 1.00 1837 1914
## y0[19] 7.30 7.40 4.73 4.49 -0.37 14.97 1.00 1523 1932
## y0[20] 37.88 37.86 4.96 4.75 29.86 46.13 1.00 1717 1999
## m1[1] 5.07 5.05 1.68 1.66 2.32 7.83 1.01 370 382
## m1[2] 9.29 9.27 1.39 1.39 7.00 11.52 1.01 384 387
## m1[3] 13.52 13.50 1.17 1.14 11.60 15.40 1.01 440 501
## m1[4] 17.74 17.70 1.04 1.02 16.11 19.42 1.01 636 790
## m1[5] 21.97 21.93 1.04 1.02 20.33 23.73 1.01 1291 1147
## m1[6] 26.19 26.16 1.18 1.10 24.31 28.15 1.00 2051 1213
## m1[7] 30.42 30.40 1.41 1.30 28.11 32.69 1.01 1996 1328
## m1[8] 34.64 34.62 1.70 1.59 31.80 37.32 1.01 1617 1328
## m1[9] 38.87 38.86 2.02 1.91 35.50 42.07 1.00 1265 1363
## m1[10] 43.09 43.05 2.36 2.28 39.26 46.88 1.00 1054 1219
## y1[1] 4.97 4.87 4.89 4.67 -2.75 13.06 1.00 1171 1582
## y1[2] 9.30 9.38 4.76 4.61 1.35 17.06 1.00 1517 1969
## y1[3] 13.54 13.49 4.79 4.46 5.83 21.18 1.00 1407 1548
## y1[4] 17.64 17.66 4.64 4.17 9.98 25.38 1.00 1709 1547
## y1[5] 22.12 22.06 4.63 4.49 14.41 29.82 1.00 1798 1891
## y1[6] 26.15 26.17 4.79 4.64 18.23 34.27 1.00 1946 1930
## y1[7] 30.46 30.49 4.98 4.68 22.06 38.43 1.00 1942 1893
## y1[8] 34.75 34.65 4.79 4.74 27.32 42.63 1.00 2056 2057
## y1[9] 38.78 38.70 5.12 4.69 30.08 47.08 1.00 1774 1528
## y1[10] 43.32 43.36 5.17 4.97 34.90 51.55 1.00 1749 1804
#explanatory variable has noise
x10=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x10=x10)
mdl=cmdstan_model('./ex9.stan')
mcmc=goMCMC(mdl,data,wrm=500,smp=1000)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 1 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 2 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 2 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 3 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 4 Iteration: 1 / 1500 [ 0%] (Warmup)
## Chain 3 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 4 Iteration: 501 / 1500 [ 33%] (Sampling)
## Chain 1 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 1 finished in 0.5 seconds.
## Chain 2 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 3 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 4 Iteration: 1500 / 1500 [100%] (Sampling)
## Chain 2 finished in 0.5 seconds.
## Chain 3 finished in 0.5 seconds.
## Chain 4 finished in 0.5 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.5 seconds.
## Total execution time: 0.7 seconds.seeMCMC(mcmc,'[m,x]')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -42.85 -44.66 11.25 9.48 -58.31 -21.40 1.05 74 155
## b0 5.36 5.35 1.82 1.81 2.32 8.29 1.00 715 936
## b1 1.87 1.87 0.18 0.18 1.57 2.17 1.01 733 1149
## s 3.26 3.34 1.37 1.44 0.95 5.42 1.03 137 108
## sx 1.50 1.52 0.83 0.94 0.25 2.86 1.05 66 162
## x0[1] 7.85 7.65 1.65 2.01 5.66 10.58 1.03 197 1456
## x0[2] 5.46 5.37 1.06 0.88 3.76 7.24 1.01 1568 1743
## x0[3] 3.78 3.99 1.42 1.48 1.30 5.64 1.02 275 763
## x0[4] 3.17 3.29 1.12 0.92 1.23 4.85 1.01 1372 1559
## x0[5] 12.12 12.11 0.99 0.78 10.52 13.78 1.01 2989 2016
## x0[6] 13.32 13.08 1.50 1.58 11.38 15.92 1.03 246 1205
## x0[7] 2.66 2.87 1.53 1.58 -0.03 4.64 1.02 255 805
## x0[8] 10.71 10.60 1.08 0.87 9.05 12.56 1.01 1466 1883
## x0[9] 20.19 20.09 1.25 0.97 18.29 22.32 1.01 1654 1534
## x0[10] 0.77 0.63 1.23 1.15 -0.97 2.95 1.01 516 741
## x0[11] 14.10 14.08 1.05 0.87 12.42 15.84 1.02 2872 1845
## x0[12] 6.66 6.84 1.40 1.54 4.26 8.51 1.02 232 1116
## x0[13] 1.33 1.16 1.37 1.45 -0.56 3.74 1.01 294 779
## x0[14] 16.99 17.05 1.14 0.90 15.13 18.81 1.01 1736 2146
## x0[15] 8.25 8.26 1.01 0.85 6.64 9.86 1.01 2594 1872
## x0[16] 5.06 5.26 1.79 2.18 1.94 7.39 1.02 197 736
## x0[17] 3.26 3.16 1.08 0.90 1.59 5.09 1.02 1266 1574
## x0[18] 4.10 4.25 1.17 1.06 2.02 5.79 1.01 541 1097
## x0[19] 1.50 1.43 1.10 0.91 -0.26 3.38 1.01 1237 1235
## x0[20] 17.96 17.73 1.32 1.08 16.16 20.33 1.01 538 1386
## m0[1] 20.06 19.85 3.07 3.86 15.53 24.83 1.03 185 1873
## m0[2] 15.61 15.64 2.01 1.86 12.30 18.87 1.00 1582 2515
## m0[3] 12.49 12.55 2.58 2.97 8.50 16.52 1.02 284 1397
## m0[4] 11.35 11.32 2.05 1.86 8.00 14.67 1.01 2435 2354
## m0[5] 28.04 27.99 1.88 1.71 25.00 31.21 1.00 4183 2241
## m0[6] 30.26 30.13 2.74 3.14 25.98 34.49 1.02 257 1732
## m0[7] 10.41 10.53 2.74 3.20 6.17 14.68 1.02 296 1806
## m0[8] 25.41 25.46 2.01 1.82 22.10 28.64 1.01 1619 1979
## m0[9] 43.08 43.13 2.36 2.12 39.13 46.84 1.00 2611 1949
## m0[10] 6.86 6.97 2.41 2.47 2.81 10.54 1.01 417 1956
## m0[11] 31.73 31.72 2.03 1.89 28.41 34.99 1.01 3918 2699
## m0[12] 17.86 17.98 2.59 3.00 13.77 21.81 1.02 243 2358
## m0[13] 7.91 8.05 2.65 2.96 3.51 11.90 1.02 273 1589
## m0[14] 37.12 37.08 2.14 1.92 33.71 40.68 1.00 2040 2190
## m0[15] 20.82 20.82 1.91 1.63 17.73 23.92 1.01 3533 2084
## m0[16] 14.90 15.10 3.27 4.10 9.90 19.68 1.02 202 1650
## m0[17] 11.51 11.59 2.03 1.87 8.26 14.71 1.01 1263 2033
## m0[18] 13.07 13.08 2.14 2.18 9.65 16.51 1.01 530 2033
## m0[19] 8.23 8.32 2.12 1.97 4.69 11.56 1.00 1819 1870
## m0[20] 38.93 39.10 2.42 2.35 34.79 42.53 1.01 794 2214
## y0[1] 20.03 20.72 4.72 4.53 11.50 26.46 1.01 449 2597
## y0[2] 15.56 15.73 4.09 3.45 8.49 22.14 1.01 3522 2430
## y0[3] 12.43 11.81 4.44 4.03 6.03 20.52 1.01 894 2754
## y0[4] 11.24 11.10 4.06 3.44 4.73 18.26 1.00 4138 3241
## y0[5] 28.06 28.02 4.03 3.46 21.46 34.73 1.01 4322 3611
## y0[6] 30.31 30.84 4.51 4.20 22.35 36.69 1.01 712 2820
## y0[7] 10.41 9.83 4.52 4.10 3.87 18.51 1.01 784 2821
## y0[8] 25.48 25.69 4.04 3.41 18.55 31.85 1.00 3597 3159
## y0[9] 43.05 43.17 4.33 3.73 35.76 50.01 1.01 3471 2935
## y0[10] 6.79 7.26 4.27 3.72 -0.79 13.15 1.01 1116 2564
## y0[11] 31.61 31.63 4.09 3.42 24.88 38.29 1.00 3978 3262
## y0[12] 17.92 17.42 4.45 4.14 11.54 25.72 1.01 785 2850
## y0[13] 7.87 8.44 4.40 4.00 0.09 14.19 1.00 942 2023
## y0[14] 37.11 36.90 4.10 3.39 30.50 44.07 1.00 3652 3179
## y0[15] 20.90 20.88 4.00 3.52 14.47 27.42 1.01 3876 3506
## y0[16] 14.89 14.31 4.88 4.97 8.17 23.65 1.01 506 1890
## y0[17] 11.50 11.77 4.00 3.47 4.51 17.84 1.00 3109 3685
## y0[18] 13.10 12.77 4.10 3.63 6.86 20.20 1.00 2272 2881
## y0[19] 8.23 8.40 4.12 3.46 1.42 14.86 1.00 3318 3206
## y0[20] 38.95 39.40 4.31 3.67 31.16 45.41 1.00 1997 3074
## m1[1] 5.29 5.28 1.83 1.82 2.24 8.22 1.00 714 917
## m1[2] 9.46 9.48 1.51 1.52 6.91 11.88 1.00 772 896
## m1[3] 13.63 13.67 1.26 1.24 11.48 15.67 1.00 917 1200
## m1[4] 17.81 17.83 1.10 1.06 15.91 19.57 1.00 1306 1868
## m1[5] 21.98 21.98 1.08 1.04 20.15 23.75 1.00 1943 2134
## m1[6] 26.15 26.16 1.21 1.16 24.13 28.15 1.00 2015 2181
## m1[7] 30.33 30.35 1.45 1.37 27.97 32.75 1.00 1638 1864
## m1[8] 34.50 34.49 1.76 1.67 31.66 37.40 1.00 1328 1564
## m1[9] 38.67 38.69 2.10 1.97 35.32 42.09 1.00 1151 1557
## m1[10] 42.85 42.88 2.46 2.33 38.87 46.83 1.00 1043 1596
## x1[1] -0.06 -0.02 1.69 1.19 -3.03 2.70 1.02 3852 3188
## x1[2] 2.21 2.20 1.72 1.17 -0.65 4.99 1.02 3935 2477
## x1[3] 4.42 4.41 1.68 1.17 1.63 7.12 1.02 4216 2514
## x1[4] 6.65 6.66 1.70 1.17 3.81 9.56 1.01 3647 2613
## x1[5] 8.87 8.87 1.70 1.15 6.15 11.72 1.02 3522 2349
## x1[6] 11.11 11.12 1.66 1.15 8.37 13.87 1.01 3853 2517
## x1[7] 13.34 13.31 1.69 1.20 10.53 16.15 1.02 4113 3224
## x1[8] 15.60 15.57 1.69 1.19 12.75 18.48 1.02 3930 2392
## x1[9] 17.76 17.78 1.80 1.23 14.76 20.65 1.02 4020 1965
## x1[10] 20.03 20.03 1.70 1.15 17.20 22.81 1.02 3728 2793
## y1[1] 5.33 5.40 3.98 3.51 -1.29 11.94 1.00 2005 2532
## y1[2] 9.50 9.50 3.86 3.44 3.27 15.98 1.00 2301 2572
## y1[3] 13.60 13.60 3.78 3.27 7.51 19.80 1.00 3272 2985
## y1[4] 17.81 17.83 3.69 3.04 11.66 23.85 1.01 3936 3320
## y1[5] 22.00 22.00 3.69 3.01 15.94 28.10 1.01 3609 2943
## y1[6] 26.14 26.14 3.74 3.20 19.87 32.17 1.01 3420 2924
## y1[7] 30.26 30.26 3.75 3.16 24.17 36.53 1.01 3338 2443
## y1[8] 34.52 34.59 3.94 3.54 28.06 40.92 1.00 2857 2569
## y1[9] 38.70 38.72 4.05 3.62 32.07 45.26 1.00 2451 2959
## y1[10] 42.90 42.88 4.26 3.79 35.92 49.97 1.00 2165 3094
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(y-smy0$median,xlab='obs.-prd.',main='residual')
density(y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
x1=mcmc$draws('x1')
smx1=summary(x1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=smx1$median,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
ex10
outlier
y~cauchy(b0+b1*x,s)ex10.stan
data{
int N;
vector[N] x;
vector[N] y;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~cauchy(b0+b1*x,s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=cauchy_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=cauchy_rng(m1[i],s);
}
}n=20
x=runif(n,0,9)
y=rnorm(n,x*2+5,1)
x[1]=3
y[1]=25
qplot(x,y)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcition
data=list(N=n,x=x,y=y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
fn(mdl,data)## Initial log joint probability = -20554.5
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 23 -32.9092 0.000287254 0.000301992 1 1 45
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -32.91
## b0 6.32
## b1 1.82
## s 3.14
## m0[1] 11.78
## m0[2] 17.34
## m0[3] 21.92
## m0[4] 9.34
## m0[5] 22.01
## m0[6] 11.87
## m0[7] 6.80
## m0[8] 7.67
## m0[9] 16.43
## m0[10] 11.21
## m0[11] 8.21
## m0[12] 15.36
## m0[13] 13.96
## m0[14] 15.47
## m0[15] 15.20
## m0[16] 14.40
## m0[17] 16.83
## m0[18] 11.48
## m0[19] 8.46
## m0[20] 14.82
## y0[1] 10.33
## y0[2] 14.17
## y0[3] 28.26
## y0[4] 11.59
## y0[5] 14.15
## y0[6] 10.50
## y0[7] 3.89
## y0[8] 7.90
## y0[9] 18.40
## y0[10] 11.94
## y0[11] 9.42
## y0[12] 18.08
## y0[13] 11.17
## y0[14] 17.45
## y0[15] 13.38
## y0[16] 14.20
## y0[17] 22.26
## y0[18] 9.26
## y0[19] 4.56
## y0[20] 8.61
## m1[1] 6.80
## m1[2] 8.49
## m1[3] 10.18
## m1[4] 11.87
## m1[5] 13.56
## m1[6] 15.25
## m1[7] 16.94
## m1[8] 18.63
## m1[9] 20.32
## m1[10] 22.01
## y1[1] -2.32
## y1[2] 7.32
## y1[3] 12.84
## y1[4] 14.20
## y1[5] 15.47
## y1[6] 10.20
## y1[7] 11.28
## y1[8] 23.59
## y1[9] 15.70
## y1[10] 20.70
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -33.33 -33.02 1.25 1.09 -35.78 -31.94 1.01 515 899
## b0 6.49 6.48 1.54 1.44 3.89 9.10 1.02 257 299
## b1 1.78 1.78 0.34 0.33 1.23 2.34 1.01 321 465
## s 3.56 3.46 0.66 0.60 2.67 4.81 1.00 1889 1428
## m0[1] 11.85 11.81 0.85 0.84 10.50 13.28 1.01 439 734
## m0[2] 17.30 17.31 1.06 1.00 15.61 19.10 1.00 1521 1291
## m0[3] 21.79 21.75 1.74 1.67 19.02 24.68 1.01 660 1045
## m0[4] 9.46 9.42 1.10 1.05 7.66 11.34 1.01 298 299
## m0[5] 21.88 21.83 1.75 1.68 19.08 24.79 1.01 655 1067
## m0[6] 11.94 11.90 0.84 0.83 10.60 13.35 1.01 451 745
## m0[7] 6.96 6.94 1.46 1.37 4.51 9.47 1.02 259 281
## m0[8] 7.82 7.80 1.33 1.25 5.65 10.11 1.01 266 288
## m0[9] 16.40 16.42 0.95 0.90 14.86 17.98 1.00 1871 1330
## m0[10] 11.29 11.25 0.89 0.89 9.87 12.80 1.01 383 594
## m0[11] 8.35 8.32 1.25 1.18 6.32 10.48 1.01 274 311
## m0[12] 15.35 15.38 0.86 0.81 13.96 16.80 1.00 1863 1316
## m0[13] 13.99 13.97 0.79 0.77 12.73 15.31 1.00 1102 1309
## m0[14] 15.46 15.49 0.87 0.82 14.07 16.94 1.00 1894 1298
## m0[15] 15.20 15.21 0.85 0.80 13.83 16.64 1.00 1804 1359
## m0[16] 14.42 14.41 0.80 0.76 13.12 15.80 1.00 1369 1328
## m0[17] 16.80 16.83 1.00 0.95 15.21 18.48 1.00 1744 1258
## m0[18] 11.55 11.51 0.87 0.87 10.17 13.02 1.01 406 717
## m0[19] 8.59 8.57 1.22 1.15 6.58 10.65 1.01 278 311
## m0[20] 14.82 14.83 0.82 0.78 13.49 16.23 1.00 1603 1353
## y0[1] 11.81 11.81 3.77 3.55 5.88 17.92 1.00 1720 1798
## y0[2] 17.34 17.25 3.84 3.54 10.98 23.62 1.00 2106 1695
## y0[3] 21.73 21.64 3.93 3.81 15.34 28.29 1.00 1794 1688
## y0[4] 9.44 9.35 3.88 3.54 3.18 16.03 1.00 1435 1686
## y0[5] 21.92 21.95 4.02 3.87 15.50 28.58 1.00 1515 1701
## y0[6] 11.85 12.06 3.80 3.49 5.45 17.73 1.00 1844 1875
## y0[7] 6.83 6.92 3.86 3.59 0.43 12.88 1.00 1173 1769
## y0[8] 7.87 7.81 3.83 3.71 1.78 14.22 1.00 1363 1801
## y0[9] 16.38 16.31 3.69 3.63 10.56 22.56 1.00 1894 1661
## y0[10] 11.32 11.24 3.84 3.77 5.03 17.68 1.00 1981 1894
## y0[11] 8.33 8.33 3.74 3.54 2.22 14.35 1.00 1246 1677
## y0[12] 15.26 15.16 3.74 3.58 9.19 21.51 1.00 1994 1863
## y0[13] 13.96 13.83 3.85 3.80 7.63 20.37 1.00 1930 1860
## y0[14] 15.50 15.60 3.70 3.51 9.14 21.33 1.00 1906 1689
## y0[15] 15.20 15.24 3.78 3.62 8.97 21.33 1.00 1788 1638
## y0[16] 14.57 14.56 3.82 3.70 8.35 20.76 1.00 2059 1842
## y0[17] 16.74 16.76 3.76 3.49 10.35 22.84 1.00 1860 1892
## y0[18] 11.56 11.52 3.68 3.56 5.53 17.73 1.00 1862 1675
## y0[19] 8.51 8.65 3.84 3.67 1.92 14.59 1.00 1152 1798
## y0[20] 14.86 14.97 3.82 3.70 8.36 21.05 1.00 1854 1941
## m1[1] 6.96 6.94 1.46 1.37 4.51 9.47 1.02 259 281
## m1[2] 8.62 8.60 1.21 1.15 6.62 10.67 1.01 279 311
## m1[3] 10.28 10.24 1.00 0.95 8.69 11.96 1.01 324 361
## m1[4] 11.93 11.90 0.84 0.83 10.60 13.35 1.01 451 745
## m1[5] 13.59 13.57 0.79 0.79 12.34 14.90 1.00 891 1208
## m1[6] 15.25 15.26 0.85 0.80 13.88 16.69 1.00 1827 1315
## m1[7] 16.91 16.93 1.01 0.96 15.28 18.61 1.00 1706 1347
## m1[8] 18.56 18.57 1.23 1.16 16.61 20.67 1.00 1069 1118
## m1[9] 20.22 20.19 1.48 1.41 17.87 22.75 1.00 791 1065
## m1[10] 21.88 21.83 1.75 1.68 19.08 24.79 1.01 655 1067
## y1[1] 7.05 7.03 3.87 3.74 0.56 13.27 1.00 1100 1416
## y1[2] 8.68 8.68 3.74 3.63 2.36 14.62 1.00 1082 1537
## y1[3] 10.31 10.27 3.68 3.51 4.51 16.42 1.00 1565 1870
## y1[4] 11.84 11.87 3.75 3.65 5.63 18.01 1.00 1869 1932
## y1[5] 13.76 13.62 3.74 3.62 7.75 19.90 1.00 1831 1929
## y1[6] 15.25 15.21 3.78 3.75 9.19 21.44 1.00 2013 1837
## y1[7] 16.87 16.99 3.77 3.62 10.72 23.24 1.00 1908 1972
## y1[8] 18.61 18.58 3.78 3.72 12.27 24.66 1.00 1682 1777
## y1[9] 20.26 20.28 3.88 3.77 13.88 26.66 1.00 1673 1858
## y1[10] 21.94 21.99 3.97 3.92 15.30 28.12 1.00 1507 1665
mdl=cmdstan_model('./ex10.stan')
fn(mdl,data)## Initial log joint probability = -103.485
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 29 -4.37618 4.66905e-05 0.0015828 0.8264 0.8264 40
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.
## variable estimate
## lp__ -4.38
## b0 5.08
## b1 1.94
## s 0.33
## m0[1] 10.91
## m0[2] 16.84
## m0[3] 21.73
## m0[4] 8.31
## m0[5] 21.82
## m0[6] 11.00
## m0[7] 5.59
## m0[8] 6.53
## m0[9] 15.86
## m0[10] 10.30
## m0[11] 7.10
## m0[12] 14.72
## m0[13] 13.24
## m0[14] 14.84
## m0[15] 14.55
## m0[16] 13.71
## m0[17] 16.30
## m0[18] 10.58
## m0[19] 7.36
## m0[20] 14.15
## y0[1] 15.03
## y0[2] 16.86
## y0[3] 22.16
## y0[4] 7.65
## y0[5] 21.56
## y0[6] 10.75
## y0[7] 5.55
## y0[8] 6.97
## y0[9] 22.29
## y0[10] 9.87
## y0[11] 0.88
## y0[12] 14.46
## y0[13] 11.04
## y0[14] 14.75
## y0[15] 13.07
## y0[16] 13.89
## y0[17] 15.22
## y0[18] 10.77
## y0[19] 7.76
## y0[20] 13.49
## m1[1] 5.59
## m1[2] 7.40
## m1[3] 9.20
## m1[4] 11.00
## m1[5] 12.81
## m1[6] 14.61
## m1[7] 16.41
## m1[8] 18.22
## m1[9] 20.02
## m1[10] 21.82
## y1[1] 7.74
## y1[2] 6.66
## y1[3] 7.86
## y1[4] -0.33
## y1[5] 395.05
## y1[6] 11.11
## y1[7] 16.77
## y1[8] 18.20
## y1[9] 19.83
## y1[10] 22.63
## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.
##
## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -7.21 -6.77 1.50 1.12 -10.16 -5.67 1.01 528 565
## b0 5.05 5.04 0.27 0.27 4.62 5.49 1.01 456 262
## b1 1.95 1.95 0.05 0.05 1.87 2.04 1.01 580 450
## s 0.39 0.37 0.13 0.12 0.21 0.64 1.01 357 149
## m0[1] 10.90 10.90 0.15 0.16 10.65 11.13 1.00 522 403
## m0[2] 16.85 16.85 0.16 0.14 16.59 17.12 1.00 1909 1332
## m0[3] 21.77 21.75 0.26 0.22 21.37 22.20 1.00 1289 1201
## m0[4] 8.29 8.29 0.20 0.20 7.98 8.61 1.01 456 288
## m0[5] 21.86 21.85 0.26 0.23 21.46 22.30 1.00 1275 1188
## m0[6] 11.00 11.00 0.15 0.16 10.76 11.23 1.00 530 421
## m0[7] 5.56 5.56 0.26 0.26 5.16 5.98 1.01 455 264
## m0[8] 6.50 6.49 0.24 0.23 6.13 6.88 1.01 453 268
## m0[9] 15.88 15.88 0.15 0.13 15.64 16.12 1.00 1907 1452
## m0[10] 10.28 10.29 0.16 0.17 10.03 10.54 1.00 492 360
## m0[11] 7.07 7.07 0.23 0.22 6.73 7.44 1.01 453 281
## m0[12] 14.73 14.74 0.14 0.12 14.50 14.95 1.00 1751 1517
## m0[13] 13.24 13.24 0.13 0.13 13.02 13.44 1.00 946 1263
## m0[14] 14.85 14.86 0.14 0.12 14.62 15.07 1.00 1781 1493
## m0[15] 14.56 14.56 0.13 0.12 14.33 14.77 1.00 1708 1517
## m0[16] 13.71 13.72 0.13 0.12 13.49 13.92 1.00 1154 1366
## m0[17] 16.32 16.32 0.15 0.13 16.06 16.57 1.00 1922 1407
## m0[18] 10.57 10.58 0.16 0.16 10.32 10.81 1.00 504 370
## m0[19] 7.34 7.34 0.22 0.22 7.00 7.70 1.01 453 280
## m0[20] 14.15 14.16 0.13 0.12 13.92 14.36 1.00 1365 1402
## y0[1] 3.06 10.89 336.09 0.55 8.26 13.44 1.00 2015 1932
## y0[2] 15.39 16.88 56.90 0.59 14.06 19.07 1.00 2108 2055
## y0[3] 21.68 21.77 5.71 0.61 19.09 24.38 1.00 2032 2014
## y0[4] 5.20 8.26 150.26 0.60 5.80 10.90 1.00 1653 1960
## y0[5] 21.95 21.87 5.55 0.67 19.56 24.53 1.00 1980 1890
## y0[6] 11.16 11.02 8.13 0.58 8.77 13.68 1.00 2051 1822
## y0[7] 5.65 5.54 9.82 0.67 3.32 7.82 1.00 1926 2056
## y0[8] 6.11 6.47 14.94 0.64 4.38 9.03 1.00 1664 1554
## y0[9] 17.21 15.86 91.80 0.62 13.07 18.24 1.00 1926 1877
## y0[10] 10.49 10.28 19.41 0.61 8.11 12.73 1.00 1891 1932
## y0[11] 6.69 7.08 9.45 0.65 3.95 9.54 1.00 1580 1879
## y0[12] 14.80 14.74 9.15 0.57 12.48 16.92 1.00 1792 1841
## y0[13] 13.52 13.26 9.30 0.59 11.16 15.30 1.00 1674 1981
## y0[14] 14.49 14.84 23.37 0.56 12.42 17.22 1.00 2157 1958
## y0[15] 14.55 14.56 14.41 0.57 12.37 17.06 1.01 1900 1742
## y0[16] 14.40 13.71 24.40 0.58 11.41 16.28 1.00 1947 1709
## y0[17] 16.44 16.31 6.84 0.59 13.93 19.19 1.00 1926 1931
## y0[18] 10.90 10.59 16.72 0.59 8.09 13.02 1.00 1767 1873
## y0[19] 7.52 7.37 13.54 0.65 4.97 10.17 1.00 1930 1927
## y0[20] 14.70 14.15 26.56 0.58 11.80 16.63 1.00 2103 1931
## m1[1] 5.56 5.56 0.26 0.26 5.16 5.98 1.01 455 264
## m1[2] 7.37 7.37 0.22 0.22 7.04 7.73 1.01 453 280
## m1[3] 9.18 9.18 0.18 0.19 8.90 9.47 1.01 468 310
## m1[4] 10.99 11.00 0.15 0.16 10.75 11.23 1.00 529 404
## m1[5] 12.81 12.81 0.13 0.13 12.58 13.01 1.00 812 1281
## m1[6] 14.62 14.62 0.13 0.12 14.39 14.83 1.00 1724 1505
## m1[7] 16.43 16.43 0.15 0.13 16.17 16.68 1.00 1922 1407
## m1[8] 18.24 18.23 0.18 0.16 17.94 18.55 1.00 1807 1449
## m1[9] 20.05 20.04 0.22 0.19 19.71 20.43 1.00 1612 1281
## m1[10] 21.86 21.85 0.26 0.23 21.46 22.30 1.00 1275 1188
## y1[1] 7.44 5.56 81.93 0.66 3.25 8.12 1.00 1466 1614
## y1[2] 7.57 7.38 19.49 0.62 4.81 9.67 1.00 1670 1612
## y1[3] 8.82 9.21 18.35 0.58 6.87 11.80 1.00 1964 2015
## y1[4] 10.86 10.98 23.35 0.60 8.40 13.35 1.00 1978 1814
## y1[5] 13.04 12.82 9.98 0.60 10.67 15.54 1.00 1812 1853
## y1[6] 13.24 14.63 50.56 0.58 12.18 17.24 1.00 2103 1931
## y1[7] 18.83 16.43 79.66 0.58 13.72 19.16 1.00 2008 1767
## y1[8] 17.83 18.24 16.77 0.62 15.46 20.63 1.00 2016 1958
## y1[9] 19.79 20.05 27.14 0.60 17.80 23.10 1.00 1886 1659
## y1[10] 21.96 21.85 5.35 0.67 19.70 24.02 1.00 1889 1971
ex11
censored
y i=1-N, y~N(m,s)
acutualy ya i=1-Na
lower censored yl i=1-Nl, y<L, P(y<L)=cdf(L-m /s)
upper censored yu i=1-Nu, y>U, P(y>U)=ccdf(U-m /s)
cdf(z) cumulative normal density function, P((-Infinity,z],z~N(0,1))
ccdf(z) complementary CDF, P([z,Infinity),z~N(0,1))
P(y | x,m,s)=P(ya i=1-Na)* P(yl i=1-Nl)* P(yu i=1-Nu)ex11.stan
data{
int N;
vector[N] x;
vector[N] y;
real L;
int Nl;
vector[Nl] xl;
real U;
int Nu;
vector[Nu] xu;
int N1;
vector[N1] x1;
}
parameters{
real b0;
real b1;
real<lower=0> s;
}
model{
y~normal(b0+b1*x,s);
for(i in 1:Nl)
target+=normal_lcdf(L | b0+b1*xl[i],s);
for(i in 1:Nu)
target+=normal_lccdf(U | b0+b1*xu[i],s);
}
generated quantities{
vector[N] m0;
vector[N] y0;
for(i in 1:N){
m0[i]=b0+b1*x[i];
y0[i]=normal_rng(m0[i],s);
}
vector[N1] m1;
vector[N1] y1;
for(i in 1:N1){
m1[i]=b0+b1*x1[i];
y1[i]=normal_rng(m1[i],s);
}
}n0=20
x=runif(n0,0,9)
y=rnorm(n0,10+3*x,4)
L=15
y[y<L]=L
nl=sum(y==L)
U=30
y[y>U]=U
nu=sum(y==U)
n=n0-nl-nu
qplot(x,y)
xy0=tibble(x=x,y=y)
xya=filter(xy0, y>L & y<U)
xyl=filter(xy0, y==L)
xyu=filter(xy0, y==U)
n1=10
x1=seq(min(x),max(x),length.out=n1) # new data fpr predcitiondata=list(N=n,x=xya$x,y=xya$y,N1=n1,x1=x1)
mdl=cmdstan_model('./ex8-0.stan')
mle=mdl$optimize(data=data)## Initial log joint probability = -6489.61
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## 16 -11.0342 0.000321468 0.000423701 1 1 37
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -11.03
## b0 15.27
## b1 1.16
## s 2.41
## m0[1] 17.23
## m0[2] 19.64
## m0[3] 22.27
## m0[4] 21.52
## m0[5] 24.14
## m0[6] 20.56
## m0[7] 20.88
## m0[8] 17.58
## y0[1] 18.43
## y0[2] 23.74
## y0[3] 18.30
## y0[4] 20.60
## y0[5] 22.56
## y0[6] 20.35
## y0[7] 22.56
## y0[8] 18.21
## m1[1] 15.96
## m1[2] 16.97
## m1[3] 17.98
## m1[4] 18.99
## m1[5] 20.00
## m1[6] 21.01
## m1[7] 22.02
## m1[8] 23.02
## m1[9] 24.03
## m1[10] 25.04
## y1[1] 13.74
## y1[2] 15.97
## y1[3] 22.45
## y1[4] 19.45
## y1[5] 19.88
## y1[6] 20.83
## y1[7] 17.74
## y1[8] 25.53
## y1[9] 25.37
## y1[10] 30.10mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 WARNING: There aren't enough warmup iterations to fit the
## Chain 1 three stages of adaptation as currently configured.
## Chain 1 Reducing each adaptation stage to 15%/75%/10% of
## Chain 1 the given number of warmup iterations:
## Chain 1 init_buffer = 15
## Chain 1 adapt_window = 75
## Chain 1 term_buffer = 10
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 WARNING: There aren't enough warmup iterations to fit the
## Chain 2 three stages of adaptation as currently configured.
## Chain 2 Reducing each adaptation stage to 15%/75%/10% of
## Chain 2 the given number of warmup iterations:
## Chain 2 init_buffer = 15
## Chain 2 adapt_window = 75
## Chain 2 term_buffer = 10
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 WARNING: There aren't enough warmup iterations to fit the
## Chain 3 three stages of adaptation as currently configured.
## Chain 3 Reducing each adaptation stage to 15%/75%/10% of
## Chain 3 the given number of warmup iterations:
## Chain 3 init_buffer = 15
## Chain 3 adapt_window = 75
## Chain 3 term_buffer = 10
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.0 seconds.
## Chain 2 finished in 0.0 seconds.
## Chain 3 finished in 0.0 seconds.
## Chain 4 finished in 0.0 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.0 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m')## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -12.17 -11.76 1.66 1.38 -15.41 -10.34 1.01 277 317
## b0 15.11 14.99 3.89 2.89 9.42 21.60 1.02 178 131
## b1 1.19 1.21 0.80 0.62 -0.10 2.39 1.02 193 180
## s 3.67 3.27 1.62 1.16 1.97 6.57 1.00 543 571
## m0[1] 17.12 17.07 2.67 2.03 13.00 21.54 1.03 184 141
## m0[2] 19.60 19.59 1.53 1.27 17.20 21.97 1.02 318 384
## m0[3] 22.31 22.33 1.85 1.55 19.36 25.17 1.00 772 704
## m0[4] 21.53 21.53 1.55 1.33 19.12 23.90 1.01 1199 1013
## m0[5] 24.22 24.22 2.86 2.33 19.61 28.56 1.01 342 353
## m0[6] 20.54 20.51 1.40 1.18 18.39 22.71 1.01 836 776
## m0[7] 20.87 20.83 1.42 1.20 18.72 23.09 1.01 1082 939
## m0[8] 17.49 17.43 2.47 1.90 13.72 21.60 1.03 187 144
## y0[1] 17.09 17.14 4.87 3.91 9.91 24.62 1.01 678 603
## y0[2] 19.56 19.55 4.42 3.55 12.74 26.71 1.00 1306 1191
## y0[3] 22.28 22.29 4.43 3.59 15.58 29.48 1.00 1693 1674
## y0[4] 21.36 21.43 4.32 3.57 14.44 27.97 1.01 1933 1675
## y0[5] 24.04 24.09 5.00 4.04 16.38 31.58 1.00 846 922
## y0[6] 20.57 20.60 4.27 3.59 13.78 27.24 1.00 1672 1481
## y0[7] 20.86 20.90 4.26 3.66 14.10 27.42 1.00 1700 1488
## y0[8] 17.65 17.59 4.85 3.81 9.83 24.82 1.00 644 582
## m1[1] 15.83 15.75 3.44 2.59 10.73 21.56 1.02 179 135
## m1[2] 16.86 16.80 2.82 2.14 12.49 21.55 1.03 182 140
## m1[3] 17.90 17.84 2.24 1.77 14.49 21.63 1.03 191 167
## m1[4] 18.93 18.92 1.75 1.45 16.30 21.67 1.02 230 233
## m1[5] 19.97 19.96 1.44 1.22 17.71 22.26 1.02 430 502
## m1[6] 21.00 20.97 1.44 1.22 18.79 23.26 1.00 1155 989
## m1[7] 22.04 22.06 1.73 1.46 19.34 24.70 1.00 937 980
## m1[8] 23.07 23.07 2.22 1.83 19.42 26.47 1.01 511 448
## m1[9] 24.11 24.11 2.79 2.28 19.64 28.36 1.01 350 354
## m1[10] 25.15 25.17 3.41 2.72 19.48 30.30 1.01 293 302
## y1[1] 15.85 15.72 5.28 4.49 7.53 24.24 1.01 369 451
## y1[2] 16.88 16.82 4.93 4.14 9.07 24.46 1.01 475 489
## y1[3] 17.83 17.75 4.78 3.79 11.15 25.30 1.01 718 1061
## y1[4] 18.82 18.79 4.32 3.60 12.17 25.60 1.00 1064 1207
## y1[5] 19.98 19.88 4.18 3.36 13.42 26.90 1.01 1319 1296
## y1[6] 21.12 21.08 4.32 3.60 14.59 28.04 1.00 1831 1595
## y1[7] 22.15 22.05 4.49 3.76 15.22 29.18 1.00 1717 1133
## y1[8] 23.23 23.13 4.72 3.77 15.78 30.46 1.00 1669 1406
## y1[9] 24.33 24.23 4.83 3.86 16.69 32.01 1.00 703 906
## y1[10] 25.17 25.11 5.30 4.24 16.94 33.27 1.01 523 592
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
data=list(N=n,x=xya$x,y=xya$y,
L=L,Nl=nl,xl=xyl$x,
U=U,Nu=nu,xu=xyu$x,
N1=n1,x1=x1)
mdl=cmdstan_model('./ex11.stan')
mle=mdl$optimize(data=data)## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Rejecting initial value:
## Log probability evaluates to log(0), i.e. negative infinity.
## Stan can't start sampling from this initial value.
## Initial log joint probability = -183.679
## Iter log prob ||dx|| ||grad|| alpha alpha0 # evals Notes
## Error evaluating model log probability: Non-finite gradient.
## Error evaluating model log probability: Non-finite gradient.
## 26 -23.6838 0.000505193 1.26264e-06 1 1 34
## Optimization terminated normally:
## Convergence detected: relative gradient magnitude is below tolerance
## Finished in 0.1 seconds.mle## variable estimate
## lp__ -23.68
## b0 5.29
## b1 3.81
## s 5.74
## m0[1] 11.70
## m0[2] 19.60
## m0[3] 28.25
## m0[4] 25.76
## m0[5] 34.35
## m0[6] 22.61
## m0[7] 23.68
## m0[8] 12.86
## y0[1] 15.34
## y0[2] 19.96
## y0[3] 25.82
## y0[4] 18.53
## y0[5] 32.44
## y0[6] 16.79
## y0[7] 11.90
## y0[8] 13.00
## m1[1] 7.56
## m1[2] 10.87
## m1[3] 14.17
## m1[4] 17.48
## m1[5] 20.78
## m1[6] 24.09
## m1[7] 27.39
## m1[8] 30.70
## m1[9] 34.00
## m1[10] 37.31
## y1[1] 7.78
## y1[2] 25.63
## y1[3] 21.57
## y1[4] 26.13
## y1[5] 17.69
## y1[6] 23.14
## y1[7] 28.77
## y1[8] 32.05
## y1[9] 27.18
## y1[10] 44.19mcmc=goMCMC(mdl,data)## Running MCMC with 4 parallel chains...
##
## Chain 1 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 1 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 1 Iteration: 600 / 600 [100%] (Sampling)
## Chain 2 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 2 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 2 Iteration: 600 / 600 [100%] (Sampling)
## Chain 3 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 3 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 3 Iteration: 600 / 600 [100%] (Sampling)
## Chain 4 WARNING: There aren't enough warmup iterations to fit the
## Chain 4 three stages of adaptation as currently configured.
## Chain 4 Reducing each adaptation stage to 15%/75%/10% of
## Chain 4 the given number of warmup iterations:
## Chain 4 init_buffer = 15
## Chain 4 adapt_window = 75
## Chain 4 term_buffer = 10
## Chain 4 Iteration: 1 / 600 [ 0%] (Warmup)
## Chain 4 Iteration: 101 / 600 [ 16%] (Sampling)
## Chain 4 Iteration: 600 / 600 [100%] (Sampling)
## Chain 1 finished in 0.1 seconds.
## Chain 2 finished in 0.1 seconds.
## Chain 3 finished in 0.1 seconds.
## Chain 4 finished in 0.1 seconds.
##
## All 4 chains finished successfully.
## Mean chain execution time: 0.1 seconds.
## Total execution time: 0.2 seconds.seeMCMC(mcmc,'m',ch=T)## variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
## lp__ -23.99 -23.62 1.54 1.39 -26.92 -22.18 1.01 301 373
## b0 -1.29 0.36 8.64 7.57 -17.08 9.47 1.02 189 164
## b1 5.29 4.87 1.93 1.69 3.03 8.92 1.02 211 223
## s 9.60 8.62 4.26 3.29 4.92 17.78 1.01 328 454
## m0[1] 7.62 8.63 5.67 5.08 -2.67 14.83 1.02 190 184
## m0[2] 18.61 18.87 2.98 2.60 13.65 22.94 1.01 339 440
## m0[3] 30.63 29.79 4.48 3.51 25.15 39.02 1.00 568 572
## m0[4] 27.17 26.57 3.60 2.83 22.59 33.80 1.01 837 711
## m0[5] 39.11 37.57 7.14 5.76 30.52 51.88 1.01 341 457
## m0[6] 22.80 22.66 2.91 2.47 18.38 27.72 1.01 949 691
## m0[7] 24.27 23.98 3.06 2.51 19.95 29.63 1.01 1037 738
## m0[8] 9.24 10.18 5.17 4.54 -0.24 15.86 1.02 191 196
## y0[1] 7.93 9.13 11.56 9.56 -12.07 24.57 1.00 860 596
## y0[2] 18.60 19.00 10.77 8.62 1.29 35.31 1.00 1926 1184
## y0[3] 30.53 29.56 11.70 9.41 14.36 50.77 1.00 1182 716
## y0[4] 27.33 26.77 11.14 9.06 10.79 46.96 1.00 1602 1122
## y0[5] 38.98 37.43 12.40 10.41 22.03 61.62 1.00 791 915
## y0[6] 23.12 22.75 11.19 8.65 5.34 41.00 1.01 1906 1233
## y0[7] 23.83 23.78 11.18 8.84 7.57 41.72 1.00 1822 1342
## y0[8] 9.06 10.23 11.72 9.75 -10.85 25.86 1.00 793 1071
## m1[1] 1.87 3.24 7.56 6.78 -11.75 11.41 1.02 188 165
## m1[2] 6.47 7.52 6.04 5.35 -4.51 14.14 1.02 189 185
## m1[3] 11.06 11.85 4.63 4.10 2.75 16.98 1.02 195 220
## m1[4] 15.66 16.20 3.47 3.00 9.62 20.25 1.01 234 345
## m1[5] 20.25 20.39 2.86 2.51 15.70 24.69 1.01 489 517
## m1[6] 24.85 24.47 3.15 2.55 20.48 30.68 1.01 1024 700
## m1[7] 29.44 28.63 4.16 3.22 24.31 37.19 1.01 642 609
## m1[8] 34.04 32.97 5.50 4.32 27.41 44.09 1.01 433 478
## m1[9] 38.63 37.16 6.98 5.58 30.23 51.09 1.01 346 472
## m1[10] 43.23 41.34 8.54 6.93 33.10 58.59 1.01 304 416
## y1[1] 1.81 3.49 12.90 10.47 -20.83 19.31 1.01 558 605
## y1[2] 6.20 7.70 12.11 9.79 -15.05 22.50 1.00 706 724
## y1[3] 11.27 11.71 10.92 9.00 -7.24 27.80 1.01 1070 830
## y1[4] 15.22 15.61 11.07 8.58 -3.31 32.12 1.00 1534 1187
## y1[5] 20.33 20.56 10.65 8.70 4.13 37.74 1.00 1673 933
## y1[6] 24.64 24.42 11.17 8.93 8.45 42.00 1.00 1860 1549
## y1[7] 29.18 28.60 11.15 9.44 12.88 48.39 1.00 1590 1062
## y1[8] 34.79 33.42 11.95 9.52 18.17 55.57 1.00 1316 1102
## y1[9] 38.22 36.74 12.71 9.75 20.87 60.77 1.00 1145 1247
## y1[10] 43.22 41.08 13.56 10.57 25.90 67.56 1.01 499 514
y0=mcmc$draws('y0')
smy0=summary(y0)
grid.arrange(
qplot(xya$y,smy0$median,xlab='obs.',ylab='prd.')+geom_abline(intercept=0,slope=1))
par(mfrow=c(1,2))
hist(xya$y-smy0$median,xlab='obs.-prd.',main='residual')
density(xya$y-smy0$median) |> plot(xlab='obs.-prd.',main='residual')
m1=mcmc$draws('m1')
smm1=summary(m1)
y1=mcmc$draws('y1')
smy1=summary(y1)
xy=tibble(x=x1,m=smm1$median,ml5=smm1$q5,mu5=smm1$q95,yl5=smy1$q5,yu5=smy1$q95)
qplot(x,y,col=I('red'))+
geom_line(aes(x=x,y=ml5),xy,col='darkgray')+
geom_line(aes(x=x,y=mu5),xy,col='darkgray')+
geom_line(aes(x=x,y=yl5),xy,col='lightgray')+
geom_line(aes(x=x,y=yu5),xy,col='lightgray')+
geom_line(aes(x=x,y=m),xy,col='black')
bimodal distribution
mixed normal distributiondata {
int<lower=0> N; // データの数
real y[N]; // 観測データ
}
parameters {
real<lower=0, upper=1> theta; // 分布1と分布2の混合比率
real mu1; // 分布1の平均
real mu2; // 分布2の平均
real<lower=0> sigma1; // 分布1の標準偏差
real<lower=0> sigma2; // 分布2の標準偏差
}
model {
for (n in 1:N) {
target += log_mix(theta,
normal_lpdf(y[n] | mu1, sigma1),
normal_lpdf(y[n] | mu2, sigma2));
}
}y0=rnorm(100,0,1)
y1=rnorm(100,-5,1)
y2=rnorm(100,5,1)
y=sample(c(y0,y1,y2),100)
density(y) |> plot()
EM algorithm
library(mclust)
rst=Mclust(y)
summary(rst)## ----------------------------------------------------
## Gaussian finite mixture model fitted by EM algorithm
## ----------------------------------------------------
##
## Mclust E (univariate, equal variance) model with 3 components:
##
## log-likelihood n df BIC ICL
## -238 100 6 -504 -505
##
## Clustering table:
## 1 2 3
## 37 24 39rst$parameters## $pro
## [1] 0.368 0.241 0.391
##
## $mean
## 1 2 3
## -4.76533 0.00268 5.18882
##
## $variance
## $variance$modelName
## [1] "E"
##
## $variance$d
## [1] 1
##
## $variance$G
## [1] 3
##
## $variance$sigmasq
## [1] 0.819
##
##
## $Vinv
## NULLplot(rst)
